
#ifndef __QUASINEWTON_MATRIX_H__

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

/* Computes the BFGS update for Hessian approximation. */
void bfgs_hess_update(const gsl_vector *p, const gsl_vector *q, 
                      gsl_vector *Bp, gsl_matrix *B);

/* Computes the BFGS update for inverse Hessian approximation. */
void bfgs_invhess_update(const gsl_vector *p, const gsl_vector *q,
                         gsl_vector *Sq, gsl_matrix *S);

/* Computes the matrix-vector product Hg, where H is 
 * the L-bfgs inverse hessian approximation g is 
 * function gradient. A history of m previous p and q 
 * vectors specified in the mxn matrices P and Q is used. 
 * The most recent iterate is stored in the kth row. */
void lbfgs_invhess_update(const gsl_vector *p, const gsl_vector *q,
                          const gsl_matrix *P, const gsl_matrix *Q,
                          int k, gsl_vector *alpha, gsl_vector *r);

/* Computes the SR1 update for Hessian approximation 
 * (this is identical to the inverse Hessian update). */
void sr1_hess_update(const gsl_vector *p, const gsl_vector *q,
                     gsl_vector *r, gsl_matrix *B);

#define __QUASINEWTON_MATRIX_H__

#endif
